log using log_supplement.smcl, replace

********************************************************************************
*
*	The Political Geography of the January 6 Insurrectionists
*   PS: Political Science & Politics
*   Replication
*
*   Robert A. Pape
*   Kyle D. Larson
*   Keven G. Ruby
*
*	Create all tables and figures in supplement
*   
*   v2024-01-08
*
********************************************************************************

clear
estimates clear 


/*

* User Installed Packages

Jann, Ben. 2023. "ESTOUT: Stata Module to Make Regression Tables." Statistical 
Software Components, February. https://ideas.repec.org//c/boc/bocode/s439301.html.

Pisati, Maurizio. 2018. "SPMAP: Stata Module to Visualize Spatial Data." 
Statistical Software Components, January. https://ideas.repec.org//c/boc/bocode/s456812.html.

Winter, Nick. 2021. "COMBOMARGINSPLOT: Stata Module to Combine the Saved Results
from Multiple Calls to Margins into One Marginsplot." Statistical Software 
Components, November. https://ideas.repec.org//c/boc/bocode/s457804.html.
*/

/*
* spmap
ssc install spmap
* combomarginsplot
ssc install combomarginsplot
* esttab
ssc install estout
*/

* Load the data
use "PS_Jan6countydata.dta", clear

 * Table A1: Descriptive Statistics  
estpost sum ///
	wpdecline mfgempdecline mfgempdecline2 pcttrump2020 trump2020vromney	///		IVs of Interest
	pctnhwhite2020 urban distdc countypop									/// 		Consistent Controls 
	if !missing(insurrectcount)

esttab using supplement.rtf , cells("count mean sd min max") replace ///
title("Descriptive Statistics")


 * Figure A1: Correlation Plot 

pwcorr insurrectcount ///
	wpdecline mfgempdecline mfgempdecline2 pcttrump2020 pctnhwhite2020 trump2020vromney ///
	urban distdc countypop if !missing(insurrectcount)

	matrix corrmatrix = r(C)

	heatplot corrmatrix, ///
		values(format(%9.2f) size(vsmall)) ///
		color(hcl diverging, intensity(.6)) ///
		xlabel(,angle(90)) ///
		legend(off) aspect(1) lower cuts(-1.05(.1)1.05) ///  
		xsize(2) ysize(2)

graph export FigureA1.png, width(3000) replace


* Alternative Variable and Model Specifications

// Other needed variables not included in main dataset 
gen votes_trump_2020b = votes_trump_2020 / 100000
label variable votes_trump_2020b "Votes for Trump (by 100,000)"

gen logpop = log(countypop)
label variable logpop "Logged County Population, 2020"

gen pop_tot_2010b = pop_tot_2010 * 100000
label variable pop_tot_2010b "County Population, 2010 (by 100,000)"
 
foreach var of varlist wpdecline mfgempdecline mfgempdecline2 unemp_15_20 pcttrump2020 trump2020vromney pctnhwhite2020 urban distdc countypop votes_trump_2020b logpop pop_tot_2010b{
	
	egen z`var' = std(`var')
}  
 
 
* Table A2 ********************************************************************
* Additional Variable Specifications 
 
local eivs "zwpdecline"
local civs "zpctnhwhite2020 urban zdistdc"

	// Primary Model, with Electoral instead of Popular Support for Trump 
	nbreg insurrectcount 									///
	`eivs' zmfgempdecline zpcttrump2020						///
	`civs' 													///
	zcountypop												///
	, robust irr	

	estimates store a2_1

	// Primary Model, with non-imputed manufacturing data  
	nbreg insurrectcount 									///
	`eivs' zmfgempdecline2 ztrump2020vromney					///
	`civs' 													///
	zcountypop												///
	, robust irr	

	estimates store a2_2
	
	// Primary Model, with unempoyment   
	nbreg insurrectcount 									///
	`eivs' zunemp_15_20 ztrump2020vromney						///
	`civs' 													///
	zcountypop												///
	, robust irr	

	estimates store a2_3
	
	// Primary Model, with unempoyment interacted with manufacturing
	nbreg insurrectcount 									///
	`eivs' c.zunemp_15_20##c.zmfgempdecline ztrump2020vromney	///
	`civs' 													///
	zcountypop												///
	, robust irr	

	estimates store a2_4
	
	
esttab a2_1 a2_2 a2_3 a2_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A2. Alternative Variable Specifications " ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Negative binomial regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")



 * Table A3 ********************************************************************
 * Using Trump Voters as the Population Control instead of County Population 

local eivs "zwpdecline zmfgempdecline ztrump2020vromney"

	// Model with White Population Decline 
	nbreg insurrectcount 	///
	zwpdecline 				///
	`civs' 					///
	zvotes_trump_2020b	 			///
	, robust irr

	estimates store a3_1

	// Model with Manufacturing Loss
	nbreg insurrectcount 	///
	zmfgempdecline 			///
	`civs' 					///
	zvotes_trump_2020b	 			///
	, robust irr

	estimates store a3_2

	// Model with Populist Support for Trump
	nbreg insurrectcount 	///
	ztrump2020vromney 		///
	`civs' 					///
	zvotes_trump_2020b	 			///
	, robust irr

	estimates store a3_3

	// Model with All Three
	nbreg insurrectcount 	///
	`eivs' 					///
	`civs' 					///
	zvotes_trump_2020b	 			///
	, robust irr	

	estimates store a3_4

// Generate Table
esttab a3_1 a3_2 a3_3 a3_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A3. Alternative Population Control - Vote for Trump" ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Negative binomial regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")

 * Table A4 ********************************************************************
 * Using Logged Population 
 
local eivs "zwpdecline zmfgempdecline ztrump2020vromney"

	// Model with White Population Decline 
	nbreg insurrectcount 	///
	zwpdecline 				///
	`civs' 					///
	zlogpop	 			///
	, robust irr

	estimates store a4_1

	// Model with Manufacturing Loss
	nbreg insurrectcount 	///
	zmfgempdecline 			///
	`civs' 					///
	zlogpop	 			///
	, robust irr

	estimates store a4_2

	// Model with Populist Support for Trump
	nbreg insurrectcount 	///
	ztrump2020vromney 		///
	`civs' 					///
	zlogpop	 			///
	, robust irr

	estimates store a4_3

	// Model with All Three
	nbreg insurrectcount 	///
	`eivs' 					///
	`civs' 					///
	zlogpop	 			///
	, robust irr	

	estimates store a4_4

// Generate Table
esttab a4_1 a4_2 a4_3 a4_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A4. Alternative Population Control - Logged Population" ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Negative binomial regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")



 * Table A5 ********************************************************************
 * Also Including County Population in 2010
  
 
	// Model with White Population Decline 
	nbreg insurrectcount 	///
	zwpdecline 				///
	`civs' 					///
	zcountypop zpop_tot_2010b	 			///
	, robust irr

	estimates store a5_1

	// Model with Manufacturing Loss
	nbreg insurrectcount 	///
	zmfgempdecline 			///
	`civs' 					///
	zcountypop zpop_tot_2010b	 			///
	, robust irr

	estimates store a5_2

	// Model with Populist Support for Trump
	nbreg insurrectcount 	///
	ztrump2020vromney 		///
	`civs' 					///
	zcountypop zpop_tot_2010b	  			///
	, robust irr

	estimates store a5_3

	// Model with All Three
	nbreg insurrectcount 	///
	`eivs' 					///
	`civs' 					///
	zcountypop zpop_tot_2010b	 			///
	, robust irr	

	estimates store a5_4
	


esttab a5_1 a5_2 a5_3 a5_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A5. Additional Control - County Population in 2010" ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Negative binomial regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")



 * Table A6 ********************************************************************
 * Logistic Regression with Binary Specification of DV
 
	// Model with White Population Decline 
	logit insurrectbinary 	///
	zwpdecline 				///
	`civs' 					///
	zcountypop	 			///
	, robust 

	estimates store a6_1

	// Model with Manufacturing Loss
	logit insurrectbinary 	///
	zmfgempdecline 			///
	`civs' 					///
	zcountypop	 			///
	, robust 

	estimates store a6_2

	// Model with Populist Support for Trump
	logit insurrectbinary 	///
	ztrump2020vromney 		///
	`civs' 					///
	zcountypop	 			///
	, robust 

	estimates store a6_3

	// Model with All Three
	logit insurrectbinary 	///
	`eivs' 					///
	`civs' 					///
	zcountypop	 			///
	, robust	

	estimates store a6_4

 
esttab a6_1 a6_2 a6_3 a6_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A6. Alternative Model Specification - Logistic with Binary DV" ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: Logistic regression with robust standard errors. Reporting odds-ratios. Standard errors in parentheses.")



 * Table A6 ********************************************************************
 * Zero-Inflated Negative Binomial 


zinb insurrectcount ///
zwpdecline `civs' zcountypop, ///
inflate(zwpdecline `civs' zcountypop) ///
robust irr 

estimates store a7_1

// Model with Manufacturing Loss
zinb insurrectcount ///
zmfgempdecline `civs' zcountypop, ///
inflate(zmfgempdecline `civs' zcountypop) ///
robust irr 

estimates store a7_2

// Model with Absolute Measure of Support for Trump
zinb insurrectcount ///
ztrump2020vromney `civs' zcountypop, ///
inflate(ztrump2020vromney `civs' zcountypop) ///
robust irr 

estimates store a7_3

// Model with All Three
zinb insurrectcount ///
`eivs' `civs' zcountypop, ///
inflate(`eivs' `civs' zcountypop) ///
robust irr 

estimates store a7_4

esttab a7_1 a7_2 a7_3 a7_4 using supplement.rtf ///
, append se compress aic bic eform constant ///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) ///
title("Table A7. Alternative Model Specification - Zero-Inflated Negative Binomial" ) ///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() ///
lines nogaps wide ///
nodepvars ///
note("Note: ZINB regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")





* Interaction Models ***********************************************************

local counter = 1 

foreach var of varlist competitive objected BidenWonBy {

	// Model with Manufacturing Interactions 
	nbreg insurrectcount 				/// 
	c.zmfgempdecline##i.`var' 	///
	`civs' 								///
	zcountypop	 						///
	, robust irr

	estimates store i`counter'
	local counter = `counter' + 1
  	
	// Model with Populist Interactions 
	nbreg insurrectcount 				/// 
	c.ztrump2020vromney##i.`var' 	///
	`civs' 								///
	zcountypop	 						///
	, robust irr

	estimates store i`counter'
	local counter = `counter' + 1
	
}
	
esttab 		i*																				///
using supplement.rtf 																			///
, replace se compress aic bic eform constant 												///
star(* 0.05 ** .01 *** .001) onecell nobase b(%9.2f) se(%9.2f) 								///
title("Table A8. Additional Interactions" ) 												///
coeflabels(																					///		 
		_cons 				"Constant"														/// Baseline Variables
		zwpdecline			"Decline in % Non-Hispanic white, 2010-2020" 					///
		zpctnhwhite2020		"Non-Hispanic White (% Pop.)"									///
		zurban				"Medium or larger Metropolitan Area"							///
		zdistdc				"Distance to D.C. (by 1,000 km)" 								///
		zcountypop			"Population (by 100,000 people)"								///
		zpcttrump2020		"% Vote for Trump, 2020"										///
		ztrump2020vromney	"Trump Populism"												///			
		1.competitive		"Contested County"												///
		1.objected			"House Rep Objected = 1"										///
		1.BidenWonBy20		"Biden Won by 20%"												///
		zBidenAdvantage		"Biden's 2020 Margin"											///
		zmfgempdecline		"Decline in Manufacturing Employment Share (1970 vs. 2020)"		///
		zunemp_15_20 		"Average Unemployment (2015-2020)"								///		
		1.competitive#c.zwpdecline 		"Contested County * White Pop Decline" 				///	Interactions
		1.objected#c.zwpdecline 		"House Rep Objected * White Pop Decline" 			///
		1.BidenWonBy20#c.zwpdecline 	"Biden Won by 20% * White Pop Decline" 				///	
		1.competitive#c.zmfgempdecline 	"Contested County * Manuf. Decline" 				///
		1.objected#c.zmfgempdecline 	"House Rep Objected * Manuf. Decline" 				///
		1.BidenWonBy20#c.zmfgempdecline "Biden Won by 20% * Manuf. Decline" 				///	
		1.competitive#c.ztrump2020vromney 	"Contested County * Trump Populism" 			///
		1.objected#c.ztrump2020vromney 	"House Rep Objected * Trump Populism" 				///
		1.BidenWonBy20#c.ztrump2020vromney "Biden Won by 20% * Trump Populism" 				///			
			) 																				///
order() 																					///
lines nogaps wide 																			///
nodepvars 																					///
note("Note: Negative binomial regression with robust standard errors. Reporting IRRs. Standard errors in parentheses.")



 * Figure 3 - IRR by Month *****************************************************
use PS_Jan6countydata.dta, clear 
 
capture frame drop getirrs

frame create getirrs ymnum counties totperps wpd pwpd lbwpd ubwpd mfl pmfl lbmfl ubmfl

foreach num of numlist 3/28 {

su ym`num', meanonly
local ymn = r(sum)
di `ymn'

	nbreg ym`num' ///
	wpdecline mfgempdecline pcttrump2020 ///
	pctnhwhite2020 ///
	urban distdc countypop ///
	, robust irr

di `ymn'

frame post getirrs  		   /// post results to the frame
        (`num')                ///  ymnum
        (e(N))                 ///  counties = counties in the model
        (`ymn')                ///  totperps = total perps in the model
        (r(table)[1,1])        ///  wpd
        (r(table)[4,1])        ///  pwpd
        (r(table)[5,1])        ///  lbwpd = lower 95% CI for wpd
        (r(table)[6,1])        ///  ubwpd = upper 95% CI for wpd
        (r(table)[1,2])        ///  mfl
        (r(table)[4,2])        ///  pmfl
        (r(table)[5,2])        ///  lbmfl = lower 95% CI for mfgemploss
        (r(table)[6,2])        //   ubmfl = upper 95% CI for mfgemploss
}

frame getirrs {
        
        capture gen lh=50 // specifies hight of the bar label
        
        list ymnum counties totperps wpd pwpd lbwpd ubwpd mfl pmfl lbmfl ubmfl
        sort ymnum
        twoway scatter wpd ymnum, msymbol(Oh) connect(1) || ///
               rcap lbwpd ubwpd ymnum, lcolor(black%25) || ///
               scatter mfl ymnum, msymbol(Th) connect(1) lpattern(solid) || ///
               rcap lbmfl ubmfl ymnum, lcolor(black%25) || ///
               scatter lh ymnum, mlabel(totperps)          ///
               mlabposition(2) mlabangle(vertical) mlabgap(0) msymbol(none) yaxis(2) || ///
               bar totperps ymnum, fcolor(none) lwidth(vthin) ///
               yaxis(2) ytitle(none axis(2)) yscale(axis(2) range(0 5000) lstyle(none)) ylabel(, axis(2) nolab notick) ///
               yscale(axis(1) range(.9 1.2)) ylabel(0.95(.05)1.2, axis(1)) ///
               ytitle("Incidence Rate Ratio") xtitle("") ///
               xscale(range(3 28)) xmtick(#20, grid) ///
               yline(1, extend lwidth(thick) lcolor(black%20) lpattern(solid)) text(1.001 24 "no effect", size(small) color(black) place(e)) ///
               scheme(plotplain) ///
               xlabel(4 "Apr 21"  16 "Apr 22" 28 "Apr 23")   ///
               legend(order( 1 3 4 6) lab(1 "White Pop Decline") lab(3 "Manufacturing Loss") lab(4 "95% CIs") lab(6 "Cumulative Charged"))
               graph export figureA2.png, width(3000) replace  
}

log close 
translate  log_supplement.smcl log_supplement.pdf 